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Abstract. We present the results of relativistic hydrodynamic simulations of the 
collision of two dense shells in a uniform external medium, as envisaged in the 
internal shock model for BL Lac jets. The non-thermal radiation produced by 
highly energetic electrons injected at the relativistic shocks is computed following 
their temporal and spatial evolution. The acceleration of electrons at the rela- 
tivistic shocks is parametrized using two different models and the corresponding 
X-ray light curves are computed. We find that the interaction time scale of the 
two shells is influenced by an interaction with the external medium. For the cho- 
sen parameter sets, the efficiency of the collision in converting dissipated kinetic 
energy into the observed X-ray radiation is of the order of one percent. 

Key words. BL Lac objects — general; Galaxies: active - quasars; X-rays: general 
— Radio sources: general. 

1. Introduction 

BL Lac objects are thought to be dominated by relativistic jets seen at small angles to 
the line of sight ( |Urry fc Padovam| l995), and their remarkably featureless radio-through- 
X-ray spectra are well fitted by inhomogeneous jet models (Brc gman et al. 1 1987). As 
the measured spectra can be reproduced by models with widely different assumptions 
the structure of the relativistic jets remains largely unknown. Only the analysis of the 
temporal variations of the emission, and combined spectral and temporal information can 
considerably constrain the jet physics. Time scales of the observed light curves are related 
to the crossing time of the emission regions which depend on wavelength and/or the time 
Send offprint requests to: PM, e-mail: pere@mpa-garching.mpg.de 
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scales of micro-physical processes like particle acceleration and radiative losses. The mea- 
sured time lags between the light curves at different energies as well as spectral changes 
during intensity variations allow to probe the microphysics of particle acceleration and 
radiation in the jet. 

Recently, several extended observation campaigns on the prominent BL Lacs PKS 
2155-304, Mrk501, and Mrk421 by ASCA and BeppoSAX, partly simultaneously with 
RXTE and TeV telescopes, have revealed that in general the X-ray spectral index and 
the peak energy correlate well with the source intensity (for a review see Pian 2002). The 
emission of the soft X-rays is generally well correlated with that of the hard X-rays and 
lags it by 3-4 ks IjTakahashi et al. 1 1996. 2000, |Zhang et al. 1 1999, IMalizia eTaFl 2000. 
IKataoka et al. 1 2000. Fos sati et al. 1 2000). However, significant lags of both signs were 
detected from several flares (|Tanihata et al. 1 2001). From XMM — Newton observations 
of PKS 21 5-304 lEdelson et al. I f2001) give, however, an upper limit to any time lags of 
|t| < 0.3 hr. They suggest that previous claims of time lags of soft X-rays with time 
scales of hours might be an artifact of the periodic interruptions of the low-Earth orbits 
of the satellites every ~ 1.6 hours. Large flares with time scales of ~ 1 day were detected 
with temporal lags of less than 1.5 hours between X-ray and TeV energies (for Mrk421 
see [Takahash i et al. 1 2000) . For all three sources the structure function and the power 
density spectrum analysis indicate a roll-over with a time scale of the order of 1 day 
or longer (Katao ka et al. 1 2001) which seems to be the time scale of the successive flare 
events. On shorter time scales only small power in the variability is found with a steep 
slope of the power density spectrum ~ /~ < - 2 ~ 3 ' ) (Tanihata 2002). 

These results were obtained from data with a relatively low signal-to-noise ratio, in- 
tegrated over wide time intervals (typically one satellite orbit). Uninterrupted data with 
high temporal and spectral resolution can now be provided by XMM — Newton . From 
an analysis of early data taken with the XMM — Newton EPIC cameras from Mrk421, 
the brightest BL Lac object at X-ray and UV wavelengths, for the first time the evolu- 
tion of intensity variations could be resolved on time scales of ~ 100 s (Bri nkmann et af~l 
2001). Temporal variations by a factor of three at highest X-ray energies were accompa- 
nied by complex spectral variations with only a small time lag of r = 265! ^ s between 
the hard and soft photons. 

In an extensive study of all currently available XMM — Newton observations of 
Mrk421 Bri nkmann et al. I C2003) find that the source exhibits a rather complex and 
irregular variability pattern - both, temporarily and spectrally. In general, an increase 
in flux is accompanied by a hardening of the spectrum as expected from a shift of the 
Synchrotron peak to higher energies. But there are exceptions and the rate of the spectral 
changes varies strongly. The shortest variability time scales appear to be of the order of 
£ ks. The lags between the hard and soft band flux are small and can be of different 
sign. 
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Correspondingly, it is hard to deduce uniquely the underlying physical parameters 
for the radiation process from the observations. For the currently favored 'shock-in-jet' 
model for the BL Lac emission (see, for example, jSpada et al. 1 2001) this implies that we 
are seeing the emission from multiple shocks which have either largely different physical 
parameters or that we detect the emission from similar shocks at very different states 
of their evolution, being additionally masked by relativistic beaming and time dilatation 
effects. 

The internal shock scenario assumes that an intermittently working central machine 
produces blobs of matter moving at different velocities along the jet. The interaction 
of two blobs is modeled as the collision of two shells whose interaction starts from the 
time of collision (Sik ora et al. I 2001; Spada et al. |Spada et al. | 2001; IModerski et al. I 
2003; Tanih ata et al. I 2003V This time is estimated from the relative velocity of two 
shells. During the interaction an internal shock propagates through the slower shell 
and accelerates electrons which produce the observed radiation ( Sp ada et al. | 2001, 
|Bicknell fc W agner 2002). These analytic models can be used to constrain the dimen- 
sions and physical properties of the emitting regions, but they cannot take into account 
the detailed hydrodynamic evolution of the interacting shells, nor the influence of the 
external medium prior to the interaction. To this end we have simulated the two dimen- 
sional axisymmetric evolution of two dense shells moving at different collinear velocities 
through a homogeneous external medium. 

In §|21 we describe the numerical method we have used to simulate both the hydro- 
dynamics and the temporal evolution of non-thermal electrons in the fluid. The shock 
acceleration process is modeled by two different approaches which are described in detail 
in §0 The results of our study are presented and discussed in §0] and the conclusions 
are given in §[S] 

2. Numerical method 

We assume that the dynamics of blazars is dominated by the thermal (baryonic or cold) 
matter while their emission is produced by a non-thermal component, which is in agree- 
ment with Siko ra et al. I (2001). This is justified in our model because the number den- 
sities of electrons and protons are equal and, thus, the inertia of the baryons is much 
larger than that of the leptons. 

We have performed a set of two dimensional axisymmetric simulations (in cylindrical 
coordinates r and z) of dense shells of matter moving at different relativistic speeds in the 
same direction. The shells collide after some time giving rise to internal shocks, where 
part of the internal energy of the thermal fluid is transferred to relativistic electrons 
producing the observed synchrotron radiation. 
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The problem is split into two parts, a thermal and non-thermal one. The evolution 
of the thermal or hydrodynamic component of rest mass density p, pressure P, radial 
velocity v r and axial velocity v z is simulated by means of a relativistic hydrodynamic 
code. The code also includes a set of N additional fluid components tracing the evolution 
of non-thermal relativistic electrons at different energies. The tracer fluids of number 
density rii (i = 1, . . . , N) are advected by the thermal fluid, and are coupled in energy 
space by their radiative energy losses. 

In the following subsections we detail the algorithms used for the simulation of the 
thermal fluid fij |2.1| L of the relativistic electrons (§ 12.211 and of the coupling between them 
(§& 



2.1. Hydrodynamics 

The equations of relativistic hydrodynamics can be cast in a system of conservation laws 
of the form 

k=l 1 

where U and F k are the vector of conserved variables and the flux vectors, respectively 
(e.g., IMarti et al. I 1994). In the case of axial symmetry and expressed in cylindrical 
coordinates (r, z) Eq. reads 

dU IdrF 8G 

at r or oz 

where 

U = [pT, phT 2 v r , phT 2 v z , phT 2 ~P-pT, pXTf (3) 
is the vector of conserved quantities and 

F = [pTv r , phT 2 v 2 + P, phT 2 v z v r , phY 2 v r - pYv r , pXTv r ] T (4) 
and 

G = [pTv Zl phT 2 v r v Zl phY 2 v 2 + P, phT 2 v z - pTv Zl pXTv z } T ', (5) 

are the corresponding flux vectors in r and z direction, respectively. Note that for reasons 
of convenience the speed of light, which is denoted by c elsewhere, is set equal to one in 
this subsection. Consequently, the Lorentz factor is given by T = (1 — v 2 — v 2 )~ x / 2 , and 
the specific enthalpy by h = 1 + e + P/p, where e is the specific internal energy of the 
fluid. 

S = [0,-,0,0,O] T (6) 

r 

is the source vector expressing the non-conservation of the radial momentum in cylindrical 
coordinates, and 



X=[X U ...,X N } T , (7) 
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is an N-component vector with 

X i = ^,i=l,...,N, (8) 
P 

where m e is the electron rest mass, and Xi is the mass fraction of the tracer species i 
with respect to the mass of the thermal fluid. 

The conservation laws are integrated with the GENESIS code of |Aloy et al. | (1999) . 
This code, which is exploits the piecewise parabolic method IjColella <k Woodw ard 1984), 
was suitably modified to passively advect a set of non-thermal species along with the main 
thermal fluid. We assume that the thermal component is a perfect fluid obeying an ideal 
equation of state of the form P = (7 ac j — l)ep, where 7 a( j = 4/3 is the adiabatic index. 

2.2. Non-Thermal population 

The non-thermal particle population evolves both in space and time. Limiting the phys- 
ical conditions in the fluid in such a way that during a time step non-thermal particles 
are contained in a single numerical cell, it is possible to split the evolution of the non- 
thermal particles in space and in time. This implies that there is a minimum magnetic 
field or, equivalently, a maximum allowed Larmor radius that depends on the numerical 
resolution employed (see below). We point out, however, that there are other approaches 
which treat the evolution of non-thermal particles by solving the diffusion-convection 
equation : Miniat i 2001. IJones et al~1 1999). The spatial evolution of the non-thermal 
particles is done by advecting them along with the fluid and, thus, their macroscopic 
velocity field corresponds to that of the thermal fluid (see § I2.1JI . The treatment of the 
temporal evolution is described in this section. 

The non-thermal particle population is assumed to be composed of ultra-relativistic 
electrons injected at shocks. Its temporal evolution is governed by the kinetic equation 
(e.g., EEEM 1962): 

^M + A ( ^ n(7 , t )) = Q (7)j (9) 

where ^(7, t) is the number density of electrons having a Lorentz factor 7 at a time t, 
and 7 = d"f/dt represents the radiative losses. In the case of synchrotron radiation the 
energy loss in cgs-units is (e.g., Ry bicki fc Lightman| 1979) 

7 = -qB 2 j 2 

with 

2e 4 



3m|c 5 

Here, B and e are the magnetic field strength and the electron charge, respectively. 

The time independent source term Q(j) present in Eq. gives the number of elec- 
trons injected at shocks with Lorentz factor 7 per unit of time. Relativistic electrons 
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are injected in zones of the computational grid which separate shocked and unshocked 
thermal fluid. Shocks are detected using the standard criterion in the Piecewise Parabolic 
Method fCo lella fc Woodw ard 1984) applied to the thermal fluid. For shock acceleration 
to take place, the magnetic field strength has to exceed some minimum value such that 
the corresponding Larmor radius tl of the fastest particle is smaller than the smallest 
zone size (AL). This is consistent with our splitting of the spatial and temporal evolution 
of the non-thermal particles (see above). Therefore, the injection of electrons at shocks 
is limited to situations where the inequality 



r L = m e c 2 V^max ~ 1 < p nm 
AL eB AL ~ [ ' 

holds. Here £ is a free parameter such that £ << 1 (in our simulations we take £ = 10 _1 ), 
and 7 max is the maximum Lorentz factor of the particles which are injected into the zone 
(see below) . In terms of the magnetic field strength, condition becomes 

^-Sat^* -1 ' (ii) 

or numerically 

1.7 • 10 3 cm\ 



B ~ { ' £AL j V^aT^lG. (12) 

The injected electrons are assumed to have a power-law distribution in the interval 
[7min7 7max] with a power law index pi ny Then the appropriate time-independent source 
term reads 

Q(l) = Qo7" P ' nJ 5( 7 ; 7mm, 7max) , (13) 
where S(x; a, b) is the interval function: 



S(x; a, b) 



1 if a < x < b 
otherwise 



The magnetic field B is assumed to be randomly oriented, and its strength is 
parametrized by the parameter as which is defined as the ratio between the energy 
density of the magnetic field and the thermal energy density of the fluid: 

B 2 p 

7T - = a B r. 

8vr 7ad - 1 

In our simulations the parameter 03 was set equal to a value of 10 -3 in order to obtain 
magnetic field strengths of the order of 0.1 G in the emitting region (Bickn eTl fc Wagner| 
2002). 

The emissivity j(u) of the synchrotron radiation for a population of relativistic elec- 
trons with a distribution 77,(7) is given by (e.g., |Rybicki fc Lightman] 1979) 

= T 2 / d7n(7)F(^/^ o7 2 ) , (14) 



(i-l)/(JV-l) 

7i =7o ( — ) , i = l,...,N 
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where F(x) = x J°° dy K 5 / 3 (y); K 5 / 3 (x) is the modified Bessel function of index 5/3, 
and 

3eB 
° 47rm c c 

In each zone of the computational grid the non-thermal electron population is repre- 
sented by a sum of power laws in N Lorentz factor intervals 

N 

n( 7 ,t) =5>' (t) 7 -««.S(7;7i-i,7i) . (15) 

i=l 

where 7i_i and ji are the lower and upper boundaries of the i-th power law distribution, 
which at time t is normalized to 71q(t;) and has a power law index Pi(t). The 7, are 
logarithmically distributed according to 

where 70 and 7at are the lower and upper limit of the whole energy interval considered 
for the non-thermal electron population. 

Eq. @ can be solved analytically for an initial power law distribution with no injec- 
tion, and for a power law injection with no initial electron distribution (Kardashcv 1962). 
We use these solutions with a slight modification arising from the fact that we solve the 
equation for each energy interval separately. 

The solution for the case of an initial power law distribution 

n(7,0) = r}o7~ Po S(7;7min,7ma X ) ( 16 ) 
of index pq with no injection (Q = 0) is 

n(7,«) = n 0l - pa (l - qB 2 1 t) p «- 2 S{ 1 /{l - qB 2 yt); 7znin, 7max) ■ (17) 

When no electrons are initially present (71(7, 0) = 0), and when the injection occurs 
at a constant rate with a power law of index p c (p c = 2.2, see § I2.4|l . i.e., 

Q(7) = Qoj- pc S( r , 7min , 7max ) (18) 

the solution is given by 

n< - 7 ' ^ = qB 2 (p° - 1) 7 ~ 2 i 1 ^" ~~ 7 high C ) 5 (^ 7min/(l + qB 2 j min t), 7 max ). (19) 
7iow and 7hi g h are evaluated according to 

(7min ) 7max ) ; if 7max/(l + gS 2 7 max t) < 7 < 7 

min 

(7min, 7/(1 - qB 2 lt)) ; if 7 nUn/(l + qB 2 j min t) < 7 < 7roin 
(7, 7/(1 - qB 2 7i) ; if 7 min < 7 < 7 max /(l + qB 2 ~f max t) 

^ (7)7max) ; if 7max/(l + qB 2 ~f max t) < 7 < 7 max 

The case 7i ow = 7, Thigh = 7/(1 — qB 2 jt) is the solution given bv IKardashevI ('1962'). 
This solution is also recovered in the limit of an infinite interval (7 m i n —> 0,7 max — > 00). 

Because the kinetic equation Eq. (jSJ) is linear in 77.(7, £)> if can b e solved for the distri- 
bution given in Eq. I(15j) using Eq. 117(1 and Eq. 1191) for each term in the sum separately. 
The new solution is then again approximated in the form of Eq. (1151) . 



(7low ) 7high ) 
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2.3. Synchrotron radiation 

The frequency dependent synchrotron emission is computed by substituting 71(7) from 
Eq. (|15fl into Eq. (|14fl . Since synchrotron self-absorption is not significant in the frequency 
range considered (10 16 -10 19 Hz), we compute the observed radiation by summing up the 
contributions from all zones of the computational grid in each time step taking into 
account the light travel time to the observer. 



2.4. Coupling of the thermal and non-thermal components 

In this section we explain how the energy losses of the non-thermal particles (see S I2.3|) 
are coupled to the thermal fluid. There are different ways in which Qq can be computed 
from the macroscopic hydrodynamic quantities. It is important to point out that the 
acceleration time scale of the electrons (jBednarz fc Ostrowskil 1996) is much smaller 
than the hydrodynamic time scale. Therefore, following the arguments of IJones et al. I 
(1999), we do not treat the shock acceleration process microscopically, but instead provide 
macroscopic models which describe the effects of the electron injection averaged over a 
hydrodynamic time step. Such models are not unique, and they depend on a number of 
free parameters, i.e., they may produce different light curves from the same hydrodynamic 
evolution. We will profit from this lack of uniqueness by comparing our synthetic light 
curves with actual observations. This will allow us to disregard injection models which 
do not match observed data. In the following two subsections we consider two models of 
electron acceleration, each with different choices of free parameters. 



2.4.1. Injection model of type-E 

Once a shock is detected, we compute e, the change of the internal energy of the 
thermal fluid per unit of time behind shocks. We assume that e aC c, the energy den- 
sity available per unit of time to accelerate non-thermal electrons, is a fraction a c of e 
flDaigne fc Mochkovitch| l998, |Bykov fc Meszaros| 1996), i.e., 

e aC c = a c e = — ^— -p , (20) 
7ad - 1 

where p denotes the temporal change of the fluid pressure behind a shock due to the 
hydrodynamic evolution. From the definition of the source term fEq. H3f) follows 

' 7maX d 77 m cC 2 Qf 7 -P- . (21) 

The back reaction of the energy loss due to particle acceleration on the flow is in- 
corporated by decreasing the pressure in the zone(s) where the acceleration takes place 
during a time interval At. ^From Eq. (|2C)|) one gets 

Ap= -2^1e acc At , (22) 
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where At is computed in the rest frame of the zone. 
Combining Eqs. l(T3j) . ({?Tjl and ({20|) one obtains 

m c^(7ad - 1)(1 - T 
where r\ = 7 max /7min- 

The three free parameters of this injection model are a e , r\ and 7 m i n , respectively. 



2.4.2. Injection model of type-N 

This model similar to the previous one in that the source term normalization Qq satisfies 
the equation l|23|) but additionally the number density of electrons n^ cc accelerated within 
a time step At is parameterized to be a fraction £ of the number of electrons in the zone 

< cc = C— , (24) 

TOp 

where m p is the proton mass. Using Eq. one finds 
(1 N 

nl cc = At— 7^7 1 f ni ( 1 - ?7 1_Pinj ) . (25) 

Pinj 

Using Eqs. J2H) and one obtains 
_ pinj - 1 a e 1 - ^-P'-j m p pAt 

7min — r, /• i 2— n ? / i\ i K^®) 

for the minimum Lorentz factor, and 

for the source normalization, where 7 m i n is given by Eq. (|26H . 

The three free parameters of this injection model are a c , and £, respectively. 



3. Results 

3.1. Hydrodynamic setup 

The density of the two colliding, identical shells is p s h = 10 4 p C xt = 10~ 22 g cm -3 , and 
their temperature is half the temperature of the external medium T cxt = 710 7 K. The 
two shells move with Lorentz factors Ti = 3 and T2 — 15, respectively. Initially, both 
shells are of cylindrical shape with a height L s h = 10 14 cm and a radius i? s h = 10 16 cm. 
The shells are initially separated by a distance Dq = 5 10 14 cm. The initial setup is shown 
schematically in Fig.^ 

The two colliding shells modeled in our simulations are initially located near the left 
boundary of the computational grid at z = z m ; n (§0. When the leading of the two shells 
approaches the right boundary of the grid at z — z max , the grid is translated into positive 
z-direction such that both shells are again located near the left boundary of the grid at 
z = z m in- In order to prevent any numerical artifact due to this re-mapping close to 



10 
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Fig. 1. A schematic view of the initial setup: Two identical, cylindrical shells of radius 
R s h, height L s h, density p s h, and initial separation D$ move through an external medium 
of density p ox t < p s h with velocities V\ and V2 along the symmetry axis (^-direction) 
such that V 2 >V\. 



the left grid boundary at z m i n , we place the shells sufficiently far from that boundary 
such that the Riemann structure emerging from the back of the trailing shell remains 
practically unaffected. 

The hydrodynamic set up consists of a two-dimensional computational grid in cylin- 
drical coordinates (r, z) of 40 x 2000 zones covering a physical domain of 1.5 10 16 cm by 
5 10 15 cm. The grid is initially filled with an external medium at rest, which has a uniform 
density p cxt = 10 -26 gcm -3 and a uniform pressure p G xt = 10~ n ergcm~ 3 . After every 
re-mapping the computational domain ahead of the shells is refilled with that external 
medium. 

We point out that the ratio \ — pc 2 /4p is rather large (according to 
|Bicknell fc Wagner] 2002) in the shells initially (x ~ 4.5 • 10 4 ), but it decreases con- 
siderably during the hydrodynamic evolution (see next section). 



3.2. Hydrodynamic evolution 

As the shells are set up as sharp discontinuities in a uniform external medium they 
experience some hydrodynamic evolution before the actual collision, which starts when 
the Riemann structure trailing the slower leading shell meets the bow shock of the fast 
trailing shell. This pre-collision evolution is rather similar for both shells, and can be 
estimated analytically using an exact one dimensional Riemann solver. In Fig.|3the two 
top panels show the analytic evolution of the flow conditions. The front (with respect to 
the direction of motion) discontinuity of each shells decays into a bow shock (5*16 and 
S2b), a contact discontinuity (in Fig.[21we only show a zoom of the one corresponding to 
the leading shell, CDIR), and a reverse shock (Sla and S2a). The back discontinuity of 
each shells develops into a rarefaction (Rib and R2b) that connects the still unperturbed 
state inside the shell with a contact discontinuity separating shell matter from the ex- 
ternal medium (in Fig. [3 only the leading shell is labelled, CD1L), and into a second 
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rarefaction that connects the contact discontinuity with the external medium (Rla and 
R2a). 

The pre-collision evolution is qualitatively similar when instead of sharp discontinu- 
ities a more smooth transition between the shells and the external medium is assumed. 
The Riemann structure emerging from the edges of the shells will be the same, i.e., it will 
consists of the same structure of shocks and rarefactions as with our set up. However, the 
exact values of the state variables in the intermediate states connecting the conditions 
in the shells with the external medium will be obviously different. 

The pre-collision hydrodynamics has two direct consequences. Firstly, each shell is 
heated by a reverse shock (Sla and S2a) which increases x by about two orders of 
magnitude (Fig.|5J). Secondly, both shells are spread in z direction as external medium 
shocked in the bow shocks (Sib and S2b) piles up in front of the shells. The latter effect 
is complicated in case of the faster trailing shell by the fact that its bow shock (5*26) 
soon starts to interact with the rarefaction (Rla) of the slower leading shell. Thereby 
the bow shock speeds up, and it eventually catches up with the slower leading shell. Our 
simulations show that the resulting interaction of the two shells occurs at a distance 
which is slightly smaller than the distance derived from an analytic estimate (see below) . 
The accelerating bow shock S2b drags along the whole Riemann structure. This explains 
why the state behind S2b is not uniform (as in case of the slower leading shell), but 
shows a monotonically decreasing density and pressure distribution (Fig.[3J). It further 
explains why the density behind the reverse shock of the faster shell (S2a) is always less 
than that behind the reverse shock (51a) of the slower shell. 

Before the bow shock S2b of the faster trailing shell can enter the interior of the slower 
shell where p — p s h, it has to cross the rarefaction Rib, i.e., it has to propagate through 
a steadily increasing density. Hence, the emission produced by the shock will increases 
gradually during this epoch until it becomes an internal shock propagating through the 
slower shell (Fig.[3] and 0J) . We point out that in the analytic model the internal shock 
does appear instantaneously when the two shells touch each other. 

Using our initial conditions (see previous section) an analytic estimate of the time 
when the two shells collide is given by 



^From our simulation we find that T w = 0.913 T c an , T 50 = 0.977 T an and T 90 = 1.009 T c an , 
where T 10 , T 50 and Tg are three moments of time at which the minimum density ahead 
of the faster shell is 10 %, 50 % and 90 % of the instantaneous maximum density p max on 
the faster shell. Note that /o max does not coincide, in general, with the initial density of 
the faster shell because the shell has undergone some hydrodynamic evolution. Although, 
pmax is only a few per cent larger than /? s h for the models under consideration. 




(28) 
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An analytic estimate of the time at which an observer located at a distance do from 
the initial position of the faster shell will receive the first light is 

T™=T™{l-V 2 /c)+t . (29) 

As mentioned in the previous section, the ratio x decreases during the interaction by 
almost four orders of magnitude because of the large pressure increase (from its initial 
value « 10 4 ; Fig. Impanel), while the density only increases by a factor of i=s 2. We 
note that the interaction region is the only one which produces a significant amount 
of all observed radiation. Hence, the value of \ m that region is the one relevant for 
observations and not the initial one. 

3.3. Setup of the non-thermal component 

The non-thermal electrons are binned into N = 48 logarithmically spaced intervals cover- 
ing a total Lorentz factor interval [71 = 1,7jv = 10 8 ] (§ 12.2(1 . The source term describing 
the electron injection at shocks (8 12.4(1 has a power law index p uli — 2.2, which is com- 
patible with the value found by Bednarz & Ostrowski (1998). The injection source terms 
are computed for the model type-E f§ l2.4.1|) with the parameters a e = 10~ 2 , i] — 10 4 and 
7min = 50, and for the model type-N f 8 12.4.2(1 with the parameters a e — 10~ 2 , rj = 10 4 
and £ = 10 -2 , respectively. 

In order to produce light curves, the synchrotron emissivity in each radiating zone is 
computed in 24 logarithmically spaced frequency bands ranging from 10 16 Hz to 10 19 Hz. 
In a subsequent step a soft and a hard light curve is obtained by integrating the com- 
puted synchrotron spectra over a prescribed energy interval. For the soft light curve we 
considered photon energies between 0.1 and IkeV, and for the hard light curve energies 
between 2 and 10 keV. 

3.4. Evolution of the non-thermal component 

FigureEl shows the normalized soft and hard X-ray light curves computed from the hy- 
drodynamic evolution using the procedure described in the previous subsection. Because 
the shells start to interact earlier than predicted by commonly used analytic estimates 
(see, e.g., |Spada et al."| 2001) the first peak of the light curve is observed well before the 
time T90 defined at the end of $ 13.21 (Fig.|5| panel b) for the injection model of type-N. 
Such an effect is not present in the injection model of type-E. The analytically predicted 
arrival time of the radiation T^ r (Eg . 129(1 should provide a good marker for the onset 
of the light curve. However, T^ r depends on the exact location of the emitting zones 
and on the specific injection model used. In the simulated collisions the spatial arrange- 
ment of the emitting zones (shocks) is different from that of non-evolving shells (analytic 



Mimica et. al: Computing X-ray light curves of blazars 13 

case), because some hydrodynamic evolution already takes place before the shells directly 
interact. The arrival time computed from our simulations is hence smaller than T™ r . 

Our simulations show that the double peak structure of the light curves is caused by 
the bow shock S2b propagating into the slower leading shell, and by the reverse shock 
S2a propagating through the faster trailing shell (see Fig.0) . Both of these two shocks are 
boosted after the actual shell interaction starts. We further find that the relative height 
of the light curve peaks depends on the injection model. In case of the model type-E 
(Fig-El panel a) the first light curve peak caused by the (former) bow shock S2b is higher 
than the second peak caused by the reverse shock S2a, because the rate of change of the 
energy density per unit of volume in shock S2b is found to be larger than in the reverse 
shock 5*20. Hence, shock S2b provides a larger amount of dissipation than shock S2a (see 
Ea. l20fl . i.e., more electrons are accelerated in the shock S2b than in the shock S2a. In 
case of the model type-N (Fig.0 panel b) we find that less electrons are accelerated in 
the shock S2b than in the case of model type-E, because their number density, which is 
proportional to the fluid density fEa. l24|) . is smaller in the (former) bow shock S2b than 
in the reverse shock 52a (see Fig©. 

Both light curves are binned into 20s time bins. Thus, if the temporal resolution 
of an observations is much worse (effectively it is ~ 80s due to signal-to-noise ratio 
limitations: IBrinkmann et al. I 2003). the finest time structures of our computed light 
curves will partially or even completely be smeared out, i.e., the differences between the 
two injection models could not be resolved. 

3.5. Energy 

Initially, the kinetic (E™ 1 ) and thermal (E^ 1 ) energies of the two shells are 

E™ = L sh R 2 sh Tr Psh (T 1 + T 2 - 2)c 2 = 4.4 1 46 erg (30) 

and 

E™ = L sh i? s 2 h 7r-^- = 4.7 1 40 erg = 9.4 lO" 7 E™ , (31) 

7ad - 1 

respectively. Beyond the time (t — 382.9 ks) at which no radiation is being produced 
anymore, the kinetic and internal energies of the fluid are 

El n = 4.35 1 46 erg = 0.99-E^ , (32) 
and 

E^ = 4.4 10 44 erg = OME^ , (33) 

respectively. The total radiated energy obtained by integrating the light curve derived 
from the integrated synchrotron spectra between 0.1 and lOkeV in time is (using model 
type-N): 

E* d = 1.4 1 43 erg- 0.03 (34) 
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with a peak luminosity of 5.8 10 39 ergs _1 . Therefore, the efficiency of converting thermal 
energy into X-ray radiation is E^ d /Ef£ ~ 0.03. 

We notice that after the collision the resulting merged shell is much hotter than the 
initial two shells (E^ >> E^ 1 ). This heating caused by the internal shocks, and by the 
pre-collision hydrodynamic evolution of the shells. 

4. Discussion 

Many of the effects found in our axisymmetric simulations can also be quite accurately 
modeled assuming a one-dimensional hydrodynamic evolution. Indeed, the initial set 
up is explicitly chosen to match this one-dimensionality, because we want to compare 
our simulation results with those obtained from previous analytic one-zone ID models. 
However, we stress that ID models represents only a first and rough attempt to model 
the physics of blazars, which is of course of multidimensional nature. In the previous 
sections we have demonstrated that our simulation results match qualitatively well with 
those of previous ID models. Thus, we are now in the position of exploring different initial 
configurations where multidimensional effects are expected to play a more important role 
(e.g., interaction of shells of different shell radii, inhomogeneities in the external medium, 
etc.). 

4.1. Shell interaction 

Concerning the temperature of the shells, which is chosen to be half the temperature of 
the external medium f$ l3.1[) . we do neither expect any significant change of the dynamic 
evolution nor the non-thermal emission, as long as density and pressure contrast between 
the shells and the external medium is sufficiently large, and as long as the temperature 
within the shells does not greatly exceed that of the ambient medium. If the latter 
requirement were to be violated, which is not expected according to observational data, 
the shells would already emit a significant amount of radiation before they collide. 

A qualitatively similar shell evolution and shell interaction is to be expected, if instead 
of two perfectly (sharp) cylindrical shells one were to consider shells having a Gaussian 
(smooth) distribution of density, pressure, and Lorentz factor. Of course, the exact val- 
ues of the arrival time, the time of interaction, etc., would be slightly different in this 
case. Decreasing the density, pressure, or Lorentz factor contrast between the shell and 
the external medium will decrease the speed of propagation of the Riemann structures 
emerging both from the leading and trailing edges (normal to the direction of propaga- 
tion) of the shell. Hence, the collision time would be even closer to the analytic estimate, 
and the shells would experience less pre-collision hydrodynamic evolution f§ l3.2[l . 

Particularly interesting would be to study (in future work) the influence of a moving 
external medium mimicking the flow of an underlying rarefied jet. In this situation the 
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Lorentz factor contrast between the shells and the external medium can be substantially 
smaller (even close to zero, if the jet moves almost as fast as the slower shell). Thus, 
the size and the depth of the Riemann structures produced by the shells will be much 
smaller. This has the important consequence that the effect of the rarefactions trailing 
the slower shell on the leading edge of the faster shell will be reduced. 

Decreasing the velocity difference between the two shells prolongs the pre-collision 
epoch, and reduces the strength of the internal shocks. Under these circumstances one 
expects that a smaller amount of energy is emitted in the X-rays bands. 

4.2. Shock acceleration models 

The time scale r acc of the shock acceleration process (r acc » 0.2s for, e.g., B = 0.1 G 
and 7 = 5 10 5 ; see § 12.4(1 is much shorter than the typical hydrodynamic time step At of 
the simulation (At w 15 s). Therefore, we have to parameterize the effect of the shock 
acceleration process on time scales of the order of the hydrodynamic time scale. 

The two different parameterizations of the shock acceleration process discussed in 
§ l2.4l produce qualitatively different light curves. An increase of the fraction of dissipated 
energy available for the acceleration of electrons (a c ) produces either more electrons 
(injection model type-E; § 12.4.1(1 . or injects electrons at higher energies (injection model 
type-N; § 12.4.2(1 . Increasing the ratio between the upper and lower Lorentz factor limit of 
the injection interval (77 = 7max/7min) produces more electrons at higher energies. Our 
particular choice of parameters is constrained by the fact that the maximum of the X-ray 
synchrotron emission is assumed to occur around 10 16 - 10 17 Hz. The particular choice of 
the parameter £ (the fraction of electrons accelerated in a zone within a hydrodynamic 
time step; 5 12.4.2(1 is suggested by fits of observed blazar spectra which imply electron 
number densities in the range of 10 3 — 10 4 cm" "3 IjMaraschi et al. 1 2003 s !. 

4.3. Efficiency 

The particular shell collision which we have simulated has an efficiency of the order of 
1 % in converting thermal energy to X-ray radiation (see 5 13.5(1 . However, if the flow 
contains more than two shells at the same time, the total efficiency of converting thermal 
to radiation energy might be larger. 

4.4. Light Curves 

The typical time scale of the flare in the observer's frame is found to be of the order of one 
kilosecond (Fig.[5j). However this finding could be influenced by our particular choice of 
the shell width and the magnetic field strength. The shape of the light curve will depend 
on the relation between the time scales associated with the radiation processes and the 
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physical size of the shells UMastichiadis fc Kirll 2002V On the one hand, the synchrotron 
cooling time scale r c of the non-thermal electrons in the comoving frame of the shell is 
t c = (qB 2 ^- 1 = 1.6 10 4 s for B = 0.25 G and 7 = 5 10 5 (see §E2J. On the other hand, 
this time scale is comparable to the light crossing time of the shell in the same frame 
(Tlc = L s h/cr ~ 3.4 10 4 s; see § I3.1|1 . Thus, electrons have sufficient time to radiate away 
a large fraction of their energy before they leave the shell. 

The light curves show a fast-rise slow-decline structure in both the soft and the hard 
energy band (Fig.lSJ). It can be recognized that in the shell collision considered in our 
simulation no significant time delays are observed between the hard and soft X-rays. 

We have checked that the patterns displayed by, at least, the hard light curves (Fig.EJ) 
are similar to the ones found bv lBrinkmann et al. 1 (2003) for Mrk421 in their Figs. 6 and 
9 (corresponding to different orbits of XMM) , respectively. We point out that the time 
scales of the observed light curves of Mrk 421 are several times longer than ours. However, 
this does not invalidate our comparison, because the duration of the synthetic light curves 
depends on the exact physical size of the shells (the wider they are in z-direction, the 
longer the duration). 

5. Conclusions 

We have performed relativistic hydrodynamic simulations of dense shells of plasma mov- 
ing in a rarefied medium, and computed the X-ray light curves produced during their 
collision. Our simulations improve existing analytic (one-zone) models (e.g., |Spada et aT~| 

2001) by computing a more realistic hydrodynamic evolution of the fluid. 

The physical conditions in the external medium correspond to the standard jet con- 
ditions in blazars IjTakahashi et al.1 2000). while those of the shells, in particular the 
density contrast between the shells and the external medium, are in agreement with the 
estimates of |Bicknell fc Wagner| (2002). These two facts enable us to properly address the 
physics of internal shocks in blazars (in contrast to what is claimed in |Bicknell fc Wagner| 

2002) . Indeed, we have demonstrated the ability of our method to produce synthetic light 
curves from relativistic hydrodynamic simulations including the back reaction of the emit- 
ted radiation on the thermal fluid. We find that both the total radiated energy and the 
light curves agree well with observational data. With our method the total efficiency of 
conversion of thermal energy into radiation energy needs not to be assumed, but can be 
directly computed. For the physical parameters used in our simulations we obtain a total 
efficiency of about 1 %. 

Our results show that the detailed structure of the synthetic light curves depends 
on the particular choice of the macroscopic parameterization of the process of particle 
acceleration. By comparing our results with observed light curves of Mrk 421 we find that 
both models of particle acceleration used in our simulations provide light curve patterns 
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which can be identified in observed light curves. However, it seems that our injection 
model type-E (where a prescribed fraction of the change of the internal energy of the 
thermal fluid per unit time is used to accelerate non-thermal electrons) fits both the 
soft and the hard light curves equally well, while the injection model type-N (similar to 
type-E, but where only a prescribed fraction of the available electrons is accelerated) fits 
only the hard part of the light curve. The different results obtained with the two models 
of injection are related to the fact that while the injection model type-E accounts for 
the strength of the shocks (the number of particles injected is proportional to the rate 
of change of energy per unit of volume behind the shock), the model type-N does not. 

The light curve produced by a single collision of two shells does not necessarily produce 
a single bumped flare, but can consist of several peaks as our simulations demonstrate. 
These peaks result from shocks which form after the actual interaction starts. One of 
these shocks propagates into the slower leading shell, while the second shock propagates 
into the faster trailing shell. This behavior is clearly different from that of one-zone 
models, which only predict a single peak from a single interaction. A comparison of our 
results with observed light curves of Mrk421 shows that some of the observed flares 
display a multiple peak structure that can be well accommodated by our model. As a 
further consequence we note that our model can explain variability in the light curves 
with less components (shells) than one-zone models. 

The results presented in this paper represent the beginning of a set of studies where we 
will consider different initial configurations in order to study the influence of such different 
initial conditions on the light curves of blazars. Among the possibilities arising we can 
immediately include more shells in the simulations, or sticking to two shell interactions 
we can vary the size of the shells. The latter variation, particularly, will allow us to 
include multidimensional effects by considering shells of different radial size. A further 
improvement of the method will be to include inverse Compton processes in the treatment 
of the radiation, which is of particular relevance in the study of blazars. 

The present approach relies on the assumption that the magnetic field is dynamically 
negligible, randomly oriented, and its energy density proportional to the thermal pressure. 
Thus, another step forward is to include dynamically important magnetic fields in the 
simulations. Finally, we also consider the possibility of performing 3D simulations which 
are better adapted to the configurations expected in nature (e.g., non-aligned shells, 
multidimensional shell trajectories, genuine 3D shell shapes, etc.). In fact both GENESIS 
and the newly coded radiative part are written (and tested) for 3D applications. However, 
due to the huge amount of computational power required for a 3D simulation, we have 
not yet performed such simulations yet. 
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Fig. 2. Snapshot illustrating the flow structure along the symmetry axis arising from the 
setup described in Fig.^just before the two shells start to interact (i = 27 ks). The lower 
panel shows the density (solid line) and pressure (dashed line) distribution measured in 
units of p oxt and /O cx tC 2 , respectively. The dash-dotted line gives the Lorentz factor of the 
fluid which is moving towards the right. The upper left (right) panel displays the exact 
solution of the one dimensional Riemann problem defined by the trailing (leading) edge 
of the right shell. Labeled are the two bow shocks Sib and S2b, the two reverse shocks 
Sla and 52a, the four rarefactions Rla, Rib, R2a and R2b, and (in the top panels only) 
the contact discontinuities CD1L and CD1R. 
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Fig. 3. Snapshots of the rest mass density (solid line; in units of 10~ 26 gem -3 ), pressure 
(dashed line; in units of 10~ 5 crgcm~ 3 ), and Lorentz factor of the fluid (dash-dotted 
line) along the symmetry axis of the computational grid at four different stages of the 
evolution. 
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Fig. 4. Contour plots of the logarithm of the rest mass density at the same moments of 
time as in Fig. [3] The coordinate values at both the r and z-axis are in units of 10 16 cm. 
The grid re-mapping procedure (see § I3.1[1 is the cause for the shift of the computational 
domain seen in panels (b) - (d). 
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Fig. 5. Normalized soft (solid line) and hard (dotted line) X-ray light curves in the 
observer's frame for injection of type-E (a) and type-N (b), respectively. The light curves 
are binned into 20s time bins. The vertical lines correspond to Tio (solid), T50 (dotted), 
Tgo (dashed), and T^" r (dash-dotted), respectively (for a definition of these times see 
$ 13. 2|) . Note that the time coordinate is renormalized to the time at which the count rate 
first exceeds of the count rate at maximum. 
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Fig. 6. Snapshot illustrating the flow structure along the symmetry axis after the two 
shells start to interact (t = 212 ks). The figure shows the density (solid line) and pressure 
(dashed line) distribution measured in units of p cxt and p C xtC 2 , respectively. The dash- 
dotted line gives the Lorentz factor of the fluid which is moving towards the right. Labeled 
are the two bow shocks Sib and S2b, the two reverse shocks Sla and 52a, the rarefactions 
Rla, Rib and R2b, and the contact discontinuity CD1R. At this time the Ricmann 
structure emerging from the leading edge of the faster shell is interacting with the trailing 
Riemann structure from the rear edge of the slower shell. This leads to the effective 
disappearance of the contact discontinuity CD1L and to the interchange of the position 
of Rla and S2b. 



